#include "solve.h"
#include "matrix.h"

double* down(Matrix a, double *p)       //解下三角形方程组
{
    p[0] = p[0] / a.datap[0];
    for (int i = 1; i < a.col; i++)
    {
        for (int k = 0; k < i; k++)
        {
            p[i] -= a.datap[i * a.col + k] * p[k];
        }
        p[i] = p[i] / a.datap[i * (a.col + 1)];
    }
    return p;
} 

double* up(Matrix a, double *p)         //解上三角形方程组（要求矩阵a行数不少于列数）
{
    p[a.col-1] = p[a.col-1] / a.datap[a.col*a.col-1];
    for (int i = a.col-2; i >=0; i--)
    {
        for (int k = a.col-1; k > i; k--)
        {
            p[i] -= a.datap[i * a.col + k] * p[k];
        }
        p[i] = p[i] / a.datap[i * (a.col + 1)];
    }
    return p;
}

double* idown(Matrix a, double *p)      //解主对角线元素均为1的下三角形方程组
{
    for (int i = 1; i < a.col; i++)
    {
        for (int k = 0; k < i; k++)
        {
            p[i] -= a.datap[i * a.col + k] * p[k];
        }
    }
    return p;
}

/*int main()
{
    double l[9] = {3, 0, 0, 1, 7, 0, 2, 5, 6};
    Matrix a(3, 3, l);
    double b[3] = {3, 8, 13};
    down(a, b);
    for (int w = 0; w < 3; w++)
    {
        std::cout << b[w] << std::endl;
    }
    a.trans();
    double c[3] = {11, 29, 18};
    up(a, c);
    for (int w = 0; w < 3; w++)
    {
        std::cout << c[w] << std::endl;
    }
    return 0;
}*/